Thorium dioxide is a robust nuclear fuel candidate for generation IV reactors because of its low generation of minor actinides, excellent radiation resistance and chemical stability.^[1–4] During the reactor operation, production of noble gases may affect the mechanical and thermal properties of nuclear fuels.^[5,6] One of these noble gases is helium, and most of the helium is generated by alpha decay during fuel burn-up. Because of its high diffusivity, helium tends to cluster and form bubbles, resulting in swelling of nuclear fuels.^[7–9] Therefore, it is necessary to investigate the clustering behavior of helium in nuclear fuels.

Considering the difficulties in nuclear fuel experiments, effective theoretical calculations should be performed. Density functional theory (DFT) is a reliable method to evaluate the point defect energies by atomic-scale calculation.^[10–13] In recent years, a number of computational works concerned with the bulk and defect properties of nuclear fuels have been performed.^[14–21] Zhang et al.^[14,15] investigated the mechanical and thermal properties of ThO₂ by first-principles calculations. Thompson et al.^[16] reported that the stability of noble gas atoms is related to the strain they caused in trap sites. Ma et al.^[17] studied the swelling of UO₂ induced by noble gases based on hybrid DFT. Brillant et al.^[18] studied the stability and solubility of fission products, including helium, using spin-polarized generalized gradient approximation with on-site Coulomb correction techniques. Yun et al.^[19] investigated the clustering behavior of He in UO₂, and found that He clusters affect the local mechanical properties of UO₂. Dabrowski et al.^[20] reported that diffusion of helium between two octahedral sites in UO₂ is along a polyline rather than a straight line.

In the present work, we investigate the clustering behavior of helium by choosing an octahedral interstitial site (OIS) and the Schottky defect (SD) as the trap sites. The volume change of ThO₂ and incorporation and solution energies are then calculated with increasing concentration of He atoms. Finally, the stability of He clusters is discussed from the perspective of their calculated binding energies.

Our calculations were performed using the density functional theory as implemented in the Vienna ab initio simulation package (VASP).^[22,23] The projected augmented wave method (PAW)^[24] and the generalized gradient approximation (GGA)^[25] were used. The exchange and correlation energies were calculated using the Perdew–Burke–Ernzerhof (PBE) functional.^[26] The wave functions were expanded in a plane-wave basis set with an energy cutoff of 500 eV. Since ThO₂ is a diamagnetic material,^[14] the spin polarization was not considered in the calculation. The results were also checked with spin polarized calculations, which showed no obvious differences. Due to no inclusion of occupied 5f states, the strong correlation effect of ThO₂ was negligible. It has been reported^[15,27] that the GGA approximation can give nearly correct energy information for ThO₂, and therefore the GGA+U method^[28,29] was not adopted in this work. The lattice constants and internal freedom of the unit cell were fully optimized until the Hellman–Feynman forces on the atoms were less than 0.01 eV/Å. The effective charge for each atom was calculated using Bader charge analysis.^[30]

In order to simulate the helium clusters incorporated in ThO₂, a 2× 2× 2 supercell containing 96 atoms was used in the calculation. The previous results^[27,28] have proven that a supercell of this size can make the energies sufficiently converged. Depending on the unit cell size and shape, a 2×2×2 Monkhorst–Pack sampling mesh^[31] of k-points was used. We implemented a k-mesh test for the He–ThO₂ system. The incorporation energy of He (discussed in the next section) calculated by 2×2×2 and 3×3×3 k-meshes has a little difference within the range of meV. This indicated that a 2×2×2 k-mesh is sufficient to avoid significant numerical errors in our calculations. All these calculations were checked using larger energy cutoffs and k-meshes; the results of total energy and Hellmann–Feynman forces were converged within 0.01 eV and 0.01 eV/Å, respectively. According to the previous work,^[32–34] the zero point energy (ZPE) of helium in oxides is small, which does not affect the numerical results. This can be seen in some previous studies of similar material systems. For instance, when helium interacts with O atoms in Al₂O₃, ZPE corrections are in a range of 10⁻²–10⁻³ eV.^[32] Therefore, calculations without ZPE correction were employed in this work.

ThO₂ crystallizes in a cubic fluorite structure (space group: Fm3m). Our calculated lattice constant is 5.617 Å, which agrees with the theoretical result (5.619 Å) reported by Zhang et al.^[14] and the experimental value (5.597 Å) reported by Olsen et al.^[35] To investigate the clustering behavior of helium in ThO₂, the following trap sites are considered: the octahedral interstitial site (OIS) and the Schottky defect (SD) as shown in Fig. 1. The SD clusters can be made by removing one thorium atom and two neighboring oxygen atoms from the lattice.^[36] Thus three configurations of the Schottky defect (SD₁, SD₂, SD₃) can be considered, differing by the distance between the two oxygen vacancies.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) (a) The cubic fluorite structure of ThO2. The gray and red spheres represent the thorium and oxygen atoms, respectively. The solid spheres represent the OIS. (b) Three configurations of the Schottky defect: (b) SD1, (c) SD2, and (d) SD3. The dashed spheres represent the vacancies.

It is known that helium atoms prefer to reside in octahedral interstitial sites in the perfect fuel matrix.^[37] He atoms spontaneously move to OISs after atomic relaxation from their initial positions in other interstitial sites. The diffusion barrier plays an important role in the clustering process of He. We calculated the migration energy of He by the nudged elastic band (NEB) method.^[38] The migration energy of He between two OISs is 3.80 eV, which is in agreement with the results reported by Da̧rowski et al.^[39] Considering the environment with a high temperature in nuclear fuels, we suggest that He tends to be mobile in ThO₂. In this work, we considered the clustering behavior of He atoms around an OIS with increasing concentration of He atoms. Firstly, we positioned two He atoms at the center of the edge between two oxygen atoms (Fig. 2), which is the midpoint between two interstitial sites. After relaxation, they both move to the center of an OIS and form a He dimer, as shown in Fig. 2(a).

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) The initial configurations of (a) two and (b) four He atoms and their configurations after atomic relaxation. The gray, red, and blue spheres represent the thorium, oxygen, and helium atoms, respectively.

The clustering behavior was investigated by increasing the concentration of He atoms in ThO₂. For four He atoms, a tetramer forms in an OIS, as shown in Fig. 2(b). For six He atoms, the He atoms do not form a hexamer, which we might except, but they form three dimers trapped in different OISs (see Fig. 3). This indicates that the OIS in ThO₂ does not have sufficient space to accommodate a large cluster, such as a hexamer. This result is different from the case in UO₂, where six He atoms tend to form a hexamer in the OIS and greatly displace the surrounding atoms.^[19] The largest displacement of the surrounding oxygen atoms in our work is 0.23 Å (see Table 1), while the displacement of oxygen atoms can be as large as 0.65 Å in UO₂. We suggest that the oxygen and thorium atoms in ThO₂ are less mobile than the uranium and oxygen atoms in UO₂.

Table 1.

Table 1.

Table 1. Incorporation energy of He in octahedral interstitial site for different k-meshes. . k-mesh Incorporation energy/eV 2×2×2 0.748 3×3×3 0.746 4×4×4 0.746

Table 1. Incorporation energy of He in octahedral interstitial site for different k-meshes. .

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) The initial configurations of six He atoms and the configuration after atomic relaxation. Three dimers form in different OISs.

When He atoms are introduced into the supercell structure, the volume of the structure changes according to the doping site and the number of He atoms. This volume change is given by

Table 2.

Table 2.

Table 2. Volume changes relative to the perfect supercell ΔV/V, displacement of the nearest-neighbouring atoms of He Δd, incorporation energies EInc, and binding energies EB with the increasing number of He atoms. . Number of He ΔV/V/% Δd/Å EInc/eV EB/eV 1 0.1 0.02 0.75 – 2 0.4 0.08 3.43 1.93 3 0.9 0.18 6.12 3.87 4 1.7 0.23 10.64 7.64 6 1.4 0.11 10.26 –

Table 2. Volume changes relative to the perfect supercell ΔV/V, displacement of the nearest-neighbouring atoms of He Δd, incorporation energies EInc, and binding energies EB with the increasing number of He atoms. .

To investigate the stability of the He clusters in ThO₂, we calculated the incorporation energies of He interstitials. The incorporation energy is defined as the energy required to incorporate a He atom in a pre-existing defect site,

The incorporation energy of He atoms is related to the strain.^[16] Owing to the large strain induced by large clusters, an OIS does not seem to be the energetically favorable trap site for He atoms. Thus, we also considered the Schottky defect (SD) as the trap site, which can provide more empty space for strain release. As shown in Fig. 4, three configurations of the Schottky defect (SD₁, SD₂, and SD₃) were considered, which differ by the distance between the two oxygen vacancies. He atoms tend to aggregate in the SDs after atomic relaxation, and even a He hexamer can exist in ThO₂ (see Fig. 4). The calculated ΔV/V and E^Inc of He clusters trapped in the SDs are listed in Table 3. As the concentration of He atoms in ThO₂ increases, both ΔV/V and E^Inc increase. A He hexamer in SD₁ has the largest incorporation energy of 3.35 eV, which also results in the largest volume change of 1.31%. As shown in Fig. 5, the incorporation energies of He atoms in the SDs are significantly smaller than those in OISs, suggesting that He clusters prefer to reside in SDs. The volume changes induced by He atoms in the SDs are also smaller than those in OISs. This indicates that these SDs are favorable for the release of the strain induced by incorporation of He clusters in ThO₂, especially for large clusters.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) The configurations of He trimmers in (a) SD1, (b) SD2, (c) SD3; and (d) a He dimer, (e) a He tetramer, and (f) a He hexamer in SD1.

Table 3.

Table 3.

Table 3. Volume changes ΔV/V (in %) and incorporation energies EInc (in eV) of He clusters trapped in SDs with the increasing number of He atoms. . Number of He SD1 SD2 SD3 ΔV/V EInc ΔV/V EInc ΔV/V EInc 1 0.39 0.01 0.42 0.03 0.30 0.20 2 0.50 0.34 0.48 0.33 0.36 0.32 3 0.56 0.57 0.55 0.62 0.40 0.66 4 0.84 1.53 0.88 1.64 0.75 1.67 6 1.31 3.35 1.16 2.67 0.88 2.41

Table 3. Volume changes ΔV/V (in %) and incorporation energies EInc (in eV) of He clusters trapped in SDs with the increasing number of He atoms. .

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) The average volume changes ΔV/V and incorporation energies EInc of He clusters trapped in SDs, with comparison to those in an OIS.

Considering the energy cost for formation of the SD, we also calculated the solution energy

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Solution energies of He clusters in SDs, with the increasing concentration of He atoms in ThO2. Three configurations of SD are considered.

Table 4.

Table 4.

Table 4. Solution energies ESol (in eV) and binding energies EB (in eV) of He clusters trapped in SDs with the increasing number of He atoms. . Number of He SD1 SD2 SD3 ESol EB ESol EB ESol EB 1 4.85 – 4.28 – 4.61 – 2 5.18 –1.16 4.58 –1.17 4.73 –1.18 3 5.41 –1.68 4.87 –1.63 5.07 –1.59 4 6.37 –1.47 5.89 –1.36 6.08 –1.33 6 8.19 –1.15 6.92 –1.83 6.82 –2.09

Table 4. Solution energies ESol (in eV) and binding energies EB (in eV) of He clusters trapped in SDs with the increasing number of He atoms. .

To access the possibility that isolated He atoms in OISs aggregate in one SD, we also calculated the binding energies (Table 4). All of the binding energies of He clusters in the SDs are negative, suggesting that the binding process is energetically favorable. For He clusters in SD₁, the He trimer has the lowest binding energy and is the most stable cluster. For SD₂ and SD₃, a He hexamer has the lowest binding energy, indicating that He atoms tend to form a large cluster in SD₂ and SD₃.

We have performed first-principles calculations of He clustering in ThO₂. As the concentration of He atoms in ThO₂ increases, the He atoms tend to form a cluster around an OIS. A He dimer is the most stable cluster in an OIS. However, one OIS is not large enough to accommodate a large cluster, such as a He hexamer. When He atoms are located in SDs, the strain induced by the He atoms is released and the incorporation and binding energies decrease. The negative binding energies indicate that He atoms located in isolated OISs can easily aggregate in a SD. A He trimer is the most stable cluster in SD₁. Large He clusters, such as a He hexamer, can also form in SDs. For SD₂ and SD₃, even large clusters (more than six He atoms) can exist according to the calculated binding energies. Finally, our results suggest that the growth of a larger He clusters may occur by the diffusion and aggregation of Schottky defects with He atoms. Our further studies may concern the formation and diffusion of these large defects in ThO₂ with a calculation using a larger supercell. The clustering behavior of He atoms will affect the mechanical properties of ThO₂. The degradation of mechanical properties will also be investigated in further investigations.